function [N, consistent, ToDo] = IA_Network_Path_Consistency_Parallel_Iteration(N, i, j, ToDo)
  consistent = 1;
  
  HasGreaterIJ = IA_Relation_Has_Greater(N(i, j));
  HasLessIJ = IA_Relation_Has_Less(N(i, j));
  HasContainsIJ = IA_Relation_Has_Contains(N(i, j));
  HasDuringIJ = IA_Relation_Has_During(N(i, j));
  
  for k = 1 : size(N, 2)
    if (k ~= i) && (k ~= j)
      % exclude equal rows
      % intervals in the same column
      if (IA_Relation_Has_Less(N(k, i)) && HasGreaterIJ)...
      || (IA_Relation_Has_Greater(N(k, i)) && HasLessIJ)...
      || (IA_Relation_Has_During(N(k, i)) && HasContainsIJ)
        continue;
      end
      % R(k,j) <- R(k,j) I N(k,j) I Constraints(N(k,i), R(i,j))
      temp = IA_Relation_Intersect(N(k, j),... 
          IA_Relation_Constraints_Hogge(N(k, i), N(i, j), N(k, j)));
      % Check consistency  
      if IA_Relation_Is_Null(temp)
        consistent = 0;
        return;
      end
      % if R(k,j) C N(k,j)
      if IA_Relation_Proper_Subset(temp, N(k, j))
        % then add <k, j> to ToDo;
        ToDo(size(ToDo, 2) + 1) = struct('i', k, 'j', j);
        N(k, j) = temp;
        N(j, k) = IA_Relation_Invert(N(k, j));
      end
    end %if
  end %for
        
  for k = 1 : size(N, 2)
    if (k ~= i) && (k ~= j)
      % intervals in the same row
      if (HasLessIJ && IA_Relation_Has_Greater(N(j, k)))...
      || (HasGreaterIJ && IA_Relation_Has_Less(N(j, k)))...
      || (HasDuringIJ && IA_Relation_Has_Contains(N(j, k)))
        continue;
      end
      % R(i,k) <- R(i,k) I N(i,k) I Constraints(R(i,j), N(j,k))
      temp = IA_Relation_Intersect(N(i, k),...
          IA_Relation_Constraints_Hogge(N(i, j), N(j, k), N(i, k)));
      % Check consistency
      if IA_Relation_Is_Null(temp)
        consistent = 0;
        return;
      end
      % if R(i,k) C N(i,k)
      if IA_Relation_Proper_Subset(temp, N(i, k))
        % then add <i, k> to ToDo;
        ToDo(size(ToDo, 2) + 1) = struct('i', i, 'j', k);
        N(i, k) = temp;
        N(k, i) = IA_Relation_Invert(N(i, k));
      end
    end %if
  end %for
end %function